function sf = solve_forward_counter(param,fp,fp_parfor,h)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Simulate the model for the counterfactuals;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

alpha = param(1:9);
TFP = param(10:11);
gamma = param(15:19);

sim_data = NaN( 3*fp_parfor.n_simulation , 13 ) ;
sim_data_network = NaN( 3*fp_parfor.n_simulation , 5 ) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Simulate the model;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

grid =  cell(fp.periods);
grid_peers =  cell(fp.periods);
grid_inv0   = cell(fp.periods);
grid_inv1   = cell(fp.periods);
hc = cell((fp.periods-1),1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
% Draw Initial Conditions;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

init_draw_log_skills = fp_parfor.init_draws;
hc{1}= exp(init_draw_log_skills);

Peers = NaN(size(hc{1},1),fp.periods);
distrib_skills   = NaN(size(hc{1},1),fp.periods);
distrib_PS = NaN(size(hc{1},1),fp.periods-1);
Skills = NaN(size(hc{1},1),fp.periods);
Peers_log = NaN(size(hc{1},1),fp.periods);
Skills(:,1) = hc{1};

%Keep track of who is moved;
moved_parent = 0.*(fp_parfor.moved_children==0) + 1.*(fp_parfor.moved_children==1);


for t=1:(fp.periods-1)

    %Random graph generation;

    [H1, H2] = ndgrid(hc{t}, hc{t});

    if t==1

        U_links1 = gamma(1).*ones(size(H1)) +  ...
            gamma(2).*log(H1) + gamma(3).*log(H2) + gamma(4).*(log(H1) - log(H2)).^(2);
        U_links2 = gamma(1).*ones(size(H1)) +  ...
            gamma(2).*log(H2) + gamma(3).*log(H1) + gamma(4).*(log(H1) - log(H2)).^(2) ;

    else

        %t-1 PS saved in the loop;
        [PS1, PS2] = ndgrid(PS, PS);
        U_links1 = gamma(1).*ones(size(H1)) +  ...
            gamma(2).*log(H1) + gamma(3).*log(H2) + gamma(4).*(log(H1) - log(H2)).^(2) + ...
            gamma(5).*( (log(H1) - log(H2)).^(2) ).*(log(H1) - log(H2)>0).*PS1;
        U_links2 = gamma(1).*ones(size(H1)) +  ...
            gamma(2).*log(H2) + gamma(3).*log(H1) + gamma(4).*(log(H1) - log(H2)).^(2) + ...
            gamma(5).*( (log(H1) - log(H2)).^(2) ).*(log(H2) - log(H1)>=0).*PS2;

    end


    Prob_Links = (1./(1+1./exp(U_links1))).*(1./(1+1./exp(U_links2))) ;

    ADJ = ( Prob_Links >= fp_parfor.peers_forward(:,:,1,t));
    ADJ(logical(eye(size(ADJ)))) = 0;

    peers_skills = hc{t} ;

    peers =  repmat(peers_skills',[size(Prob_Links,1) ,1]) ;
    temp = ADJ.*peers;
    H_peers =  nansum(temp,2)./sum(temp~=0 &  isnan(temp)==0,2) ;
    Peers(:,t) = H_peers;

    %Compute the mean-log peers for the indirect inference moments;
    temp = ADJ.*log(peers);
    H_peers_logs =  nansum(temp,2)./sum(temp~=0 &  isnan(temp)==0,2) ;
    Peers_log(:,t) = H_peers_logs;

    %Number of Friends
    n_friends = sum( ADJ , 2 );

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
    % Simulate Parental Optimal Behavior;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

    if fp.do_counter_type<3

        %type1: parents who are not moved (receiving parents);
        %type2: parents who are moved into the new neighborhood (but keep their preferences for parenting);

        State_vars = [hc{t} ,  H_peers ];
        value0_type1 = Interp_function_simulation(fp_parfor.estim_value0{t,1},State_vars);
        value1_type1 = Interp_function_simulation(fp_parfor.estim_value1{t,1},State_vars);

        value0_type2 = Interp_function_simulation(fp_parfor.estim_value0{t,2},State_vars);
        value1_type2 = Interp_function_simulation(fp_parfor.estim_value1{t,2},State_vars);

        %Probability of authoritarian parenting (by not moved (type 1) or moved (type 2);
        PS_prob_type1 = (1./(1+1./exp(value1_type1-value0_type1)));
        PS_prob_type2 = (1./(1+1./exp(value1_type2-value0_type2)));

        PS_type1 = NaN(size(value1_type1));
        PS_type2 = NaN(size(value1_type2));

        %Simulate the behavior given random numbers ( fp_parfor.PS_forward );
        PS_type1(isnan(PS_prob_type1)==0) = ( PS_prob_type1(isnan(PS_prob_type1)==0)>= fp_parfor.PS_forward(isnan(PS_prob_type1)==0,1,t) ) ;
        PS_type2(isnan(PS_prob_type2)==0) = ( PS_prob_type2(isnan(PS_prob_type2)==0)>= fp_parfor.PS_forward(isnan(PS_prob_type2)==0,1,t) ) ; 
        %Aggregate distribution of parenting style (aggregating over types);
        PS = (moved_parent==0).*PS_type1 + (moved_parent==1).*PS_type2;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
        %Optimal parental investments by parenting style and type of parent (moved vs receiving);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

        %P=0 & receiving (type1)
        inv0{1}  = Interp_function_simulation(fp_parfor.estim_policy0{t,1},State_vars);
        inv0{1}(inv0{1}<fp.Min_Time)=fp.Min_Time;
        inv0{1}(inv0{1}>fp.Max_Time)=fp.Max_Time;

        %P=0 & receiving (type2)
        inv0{2}  = Interp_function_simulation(fp_parfor.estim_policy0{t,2},State_vars);
        inv0{2}(inv0{2}<fp.Min_Time)=fp.Min_Time;
        inv0{2}(inv0{2}>fp.Max_Time)=fp.Max_Time;


        %P=1 & receiving (type1)
        inv1{1}  = Interp_function_simulation(fp_parfor.estim_policy1{t,1},State_vars);
        inv1{1}(inv1{1}<fp.Min_Time)=fp.Min_Time;
        inv1{1}(inv1{1}>fp.Max_Time)=fp.Max_Time;
        
        %P=1 & receiving (type2)
        inv1{2}  = Interp_function_simulation(fp_parfor.estim_policy1{t,2},State_vars);
        inv1{2}(inv1{2}<fp.Min_Time)=fp.Min_Time;
        inv1{2}(inv1{2}>fp.Max_Time)=fp.Max_Time;

        %Aggregate distribution of parenting style (aggregating over parenting styles and types);
        Inv = (moved_parent==0).*(inv0{1}.*(PS_type1==0) + inv1{1}.*(PS_type1==1)) + (moved_parent==1).*(inv0{2}.*(PS_type2==0) + inv1{2}.*(PS_type2==1));

        %Evolution of a child's skills given optimal parental behavior;
        Skills(:,t+1) = technology(t,alpha, hc{t}, Inv, PS, H_peers,  TFP );
        %Store distribution of skills;
        hc{t+1} =  Skills(:,t+1);

        %Updating the grid points of the model when solving backward;
        grid{t} = prctile( hc{t} , fp.percentiles_kid ) ;
        grid_peers{t} = prctile( Peers(:,t) , fp.percentiles_kid ) ;
        grid_inv0{t}  = prctile( [inv0{1} ; inv0{2} ] , fp.percentiles_inv ) ;
        grid_inv1{t}  = prctile( [inv1{1} ; inv1{2} ] , fp.percentiles_inv ) ;

        %Store the aggregate state variables (distribution of skills and
        %parenting style in the economy);
        distrib_skills(:,t+1) =  Skills(:,t+1) ;
        distrib_PS(:,t) = PS;


    elseif  fp.do_counter_type == 3 %Fixed parental behavior at baseline;

        PS  = fp_parfor.fixed_PS(:,t);
        Inv = fp_parfor.fixed_inv(:,t);
        %Evolution of a child's skills given optimal parental behavior;
        Skills(:,t+1) = technology(t,alpha, hc{t}, Inv, PS, H_peers,  TFP );
        %Store distribution of skills;
        hc{t+1} =  Skills(:,t+1);

    end


    %%%%%%%%%%%%%%%%%%%%%%%%%;
    %Store Simulated Data;
    %%%%%%%%%%%%%%%%%%%%%%%%%;
    ind_obs = 1+(t-1)*fp_parfor.n_simulation:t*fp_parfor.n_simulation;

    sim_data( ind_obs , 1 ) = 1:fp_parfor.n_simulation;
    sim_data( ind_obs , 2 ) = repmat(fp.ages(t),length(ind_obs),1);
    sim_data( ind_obs , 3 ) = h;
    sim_data( ind_obs , 4 ) = hc{t};
    sim_data( ind_obs , 5 ) = hc{t+1};
    sim_data( ind_obs , 6 ) = Inv;
    sim_data( ind_obs , 7)  = Peers(:,t) ;
    sim_data( ind_obs , 8 ) = Peers_log(:,t);
    sim_data( ind_obs , 9 ) = PS;



    sim_data_network( ind_obs , 1 ) = 1:fp_parfor.n_simulation;
    sim_data_network( ind_obs , 2 ) = repmat(fp.ages(t),length(ind_obs),1);
    sim_data_network( ind_obs , 3 ) = h;
    sim_data_network( ind_obs , 4 ) = hc{t} ;
    sim_data_network( ind_obs , 5 ) = n_friends;

end %t

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Add next-period peers into the simulated dataset;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
for t=1:(fp.periods-2)

    ind_obs = 1+(t-1)*fp_parfor.n_simulation:t*fp_parfor.n_simulation;
    sim_data( ind_obs , 10 ) = Peers(:,t+1);
    sim_data( ind_obs , 11 ) = Peers_log(:,t+1);

end

sim_data( : , 12 ) =  fp_parfor.j ;


for t=1:(fp.periods-1)

    ind_obs = 1+(t-1)*fp_parfor.n_simulation:t*fp_parfor.n_simulation;

    sim_data( ind_obs , 13 )       = fp_parfor.moved_children ;

end



sf.distrib_skills = distrib_skills;
sf.distrib_PS = distrib_PS;
sf.grid = grid;
sf.grid_peers = grid_peers;
sf.grid_inv0 = grid_inv0;
sf.grid_inv1 = grid_inv1;
sf.sim_data = sim_data;
sf.sim_data_network = sim_data_network;

